A PRESSURE BASED METHOD FOR THE SOLUTION 
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ABSTRACT , /) , 

A new technique has been developed for the solution of the incompressible 
Navier-Stokes equations. The numerical technique, derived from a pressure 
substitution method (PSM) , overcomes many of the deficiencies of the pressure 
correction method. This technique allows for the direct solution of the actual 
pressure in the form of a Poisson equation which is derived from the pressure 
weighted substitution of the full momentum equations into the continuity 
equation. Two-dimensional internal flows are computed with this method and the 
prediction of cascade performance is presented. The extension of the pressure 
correction method for the solution of three-dimensional flows is also 
presented . 
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n north face of control volume 

N north grid point 

o outlet flow conditions 

P first grid point away from the wall 

s south face of control volume 

S south grid point 

u upstream face of control volume 

U upstream grid point 


INTRODUCTION 


Turbomachinery flows are amongst the most complex encountered in 
engineering. These flows are three dimensional and highly turbulent and may 
include heat transfer. Unsteadiness and compressibility in turbomachinery 
flows are also relevant but will not be considered in this paper. Higher 
order turbulence models are also needed to resolve the transitional and near 
wall flow phenomena, which are important for loss and heat transfer 
predictions. Flow separation should also be accurately predicted with such 
models as these regions are inevitably present in modern highly loaded turbine 
and compressor blade rows. 

This paper is mainly concerned with the development of a more efficient 
pressure based method and its application to turbomachinery incompressible 
flows . 

A brief review of the relevant works follows. Hah (ref. 1) solved the 
uncoupled equations on a staggered grid. Higher order discretization was used 
for the convection terms to reduce numerical diffusion. The quadratic upsteam 
(QUICK) differencing by Leonard (ref. 2) and skew upwinding scheme by Raithby 
(ref. 3) were extended and modified. An algebraic Reynolds stress model 
modified for the effect of streamline curvature and rotation was used for the 
closure of the governing equations. Comparison with experimental data for an 
end-wall turbine cascade flow, showed that various complex three-dimensional 
viscous flow phenomena were well predicted. 

An excellent review paper by Patankar (ref. 4) puts the treatment of 
convection and diffusion into perpective. He refers to a study done by Patel 
and Markatos (ref. 5) who evaluated eight different discretization schemes to 
solve the two-dimensional convection-diffusion equations. "They found that 
QUICK and its variants gave accurate and cost effective solutions when they 
converged. However they failed to converge for high flow rates and coarse 
grids." The situation becomes more acute for supersonic flows for which Patel 
et al. (ref. 6) concluded that the upwind scheme is the best choice. 

Patankar (ref. 4) observes, "Thus, we have gone around full circle. 
Lower-order schemes such as upwind are stable and monotonic but lead to false 
diffusion. Higher-order schemes such as QUICK eliminate false diffusion but 
produce wiggles and often fail to converge. The search for the perfect method 
is still not over." 

Rhie and Chow (ref. 7) presented a solution of the Navi er-Stokes equations 
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in general curvilinear coordinates on a non-staggered grid. A specific 
algorithm based on pressure weighting (PWM) was developed to suppress the 
pressure oscillations. This development was a major breakthrough for the 
solution of the fluid flow equations in generalised coordinates. They used the 
pressure correction method (PCM) to couple the momentum and continuity 
equations. The standard form of the k-f model was used to describe the 
turbulence for flows over airfoils with and without trailing edge separation. 

Recently, Rhie (ref. 8), extended his scheme to employ the pressure 
implicit split operator (PISO) concept with a multigrid procedure to enhance 
convergence. Within the several levels of pressure correction, the density was 
treated implicitly for the correction of the mass flow imbalance. This 
treatment was applicable for compressible flows. The method was tested for a 
wide range of flows, including three-dimensional driven cavity and a turbine 
end-wall cascade. 

Three different formulations for non-staggered grids have been suggested 
by Shih and Ren (ref. 9). Some employ the Poisson equation for pressure, in 
place of the continuity equation. All of the algorithms, proposed by them, 
solve for the pressure directly and not the pressure correction as in the 
SIMPLE algorithm. When the continuity equation is used, a Pressure 
Substitution Method (PSM) is employed. However the equations were solved in 
non-conservative form, thus resulting in a purely finite difference method. 

Hobson and Lakshminarayana (ref. 10) presented a control volume 
formulation of the pressure substitution method for the solution of the 
Navier-Stokes equations describing incompressible flows. The application of 
this method to two-dimensional cascade flows is presented. The Low-Reynolds 
number k-e turbulence model according to Lam and Bremhorst (ref. 11) is used 
to resolve the turbulence within the flow. 

Three-dimensional flows are calculated with the pressure correction method 
similar to that used by Rhie (ref. 12) and the standard High-Reynolds number 
k-e turbulence model including wall functions is used. The computation of 
turbulent flow in a turbine end-wall cascade is presented. 

THEORETICAL FORMULATION 


The Incompressible Viscous Flow Equations 


For steady-mean incompressible flow, the time-averaged continuity equation 
and the Reynolds averaged Navier-Stokes equation can be written in Cartesian 
tensor form: 


f>i> = 0 
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where p is the mean density, the mean velocity, p 
the molecular viscosity. The instantaneous velocity, 
implies a time averaged quantity. 


the mean pressure and p 

u! , with an overbar " — " 

1 ’ 


For the prediction of heat transfer in non-i sothermal flows through 
turbines, the incompressible energy equation is also solved. The form used is 
given as follows, 


d(pu.h) 

dx. 1 
J 


_ d_ / r T«9T 

dx . 1 dx . 
J l i 



+ p(§_ + 0’ ) 


[3] 


T 

where h is the static enthalpy, T the thermal conductivity of the gas, T the 
static temperature. Once again, the averaged fluctuating correlations are 
present, which will be described in the next section. The dissipation 

function, <|), is 



du.ldu. 
+ Sx 1 dx-! 


J i 


[4] 


Assuming that the specific heat, C , is constant and the gas is ideal, the 

static enthalpy is replaced by the static temperature, according to the 
following relationship. 


h = 


C T 
P 


[5] 


Two-Equation Low-Reynolds Number Turbulence Model 


To solve equations [2] and [3], one must specify the Reynolds shear 

stress, pu!u! , and heat flux, pu!T’ . A turbulent or "eddy" viscosity, p, , and 
i j i t 

a turbulent Prandtl number, Pr^, are defined according to the Boussinesq 
approximation (ref. 13) such that, 


-pu!ul = p. 
i J 


"du. du.l 2, c 
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-pu:r = Pr 


M t r^T 
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dx. 


[ 7 ] 


where 6.. is the Kronecker delta and k is the turbulent kinetic energy. 

ij 

The isotropic eddy viscosity model was used to obtain an effective 
viscosity, which is the sum of the turbulent viscosity and molecular 
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viscosity. 


For the purpose of this paper, the turbulent Prandtl number will be 
assumed constant and equal to 0.9. Although this is not generally true, it has 
been found to be a reasonably good approximation for most boundary layer flows 
with mild temperature gradients. The role of the turbulence model now becomes 
that of determining the correct value of p^. 

From the k-e turbulence model, as proposed by Launder and Spalding (ref. 
14), the turbulent eddy viscosity p^ is given by, 


p. = f C p - 
't p p r e 


where k is the turbulent kinetic energy and £ is the turbulent energy 
dissipation and each are defined as follows, 


[ 8 ] 


k = 


77 u:u. 

2 l l 


[9] 


p d u! chi! 

e ~ p dx 1 . dx ) 
J J 


[ 10 ] 


The modelled equations for the transport of k and e are given as follows, 
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Upon consideration of equations [8] 
are five empirical constants; C , C , 

functions; f , f , and f . The values, 
’ p' 1 ’ 2 ’ 

paper are those suggested as "standard” 

These assume the following values, 


, [11] and [12], one can see that there 

C , a. , o , and three empirical 
£2’ k’ £' 

of the five constants, used in this 
by Launder and Spalding (ref. 14). 


C =0.09, C =1.44, C =1.92, 07=1.0, a =1 . 3 
p ’ £1 ’ £2 ’ k £ 


[13] 


Turbulent motions immediately adjacent to the wall are significantly 
influenced by the presence of the wall. Here the magnitude of the effective 
turbulent viscosity becomes small, and the effects of the molecular viscosity 
become important. Experimental work has shown that in some turbulent flow 
situations, there exists a common structure or behavior near the wall. Under 
these conditions both the mean velocities and the measurable turbulence 
quantities exhibit nearly universal behavior. The knowledge of this structure 
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has allowed the formulation and use of the so-called wall-functions. These 
functions algebraically bridge the near wall region and eliminate the need for 
more expensive and time consuming calculations with a fine grid near the wall. 

Unfortunately, there are many flow situations of interest where this 
near-wall similarity breaks down. Large pressure gradients are found in the 
leading edge and trailing edge regions of turbomachinery blade profiles due to 
the high streamline curvature in these areas. This, and mass transfer at the 
wall, for example, both alter the near-wall flow, thus wall functions cannot 
be used. To incorporate these effects into the turbulence models, a variety of 
suggestions have been made. The purpose of the functions f , f and f is to 

fJ/ 1 2 

provide a similar modification to the k-e model, thus extending the validity 
of the equations through the viscous sub-layer to the wall. Use is made of two 
"turbulent Reynolds numbers." These are defined as follows, 


Re = & 

[14] 

t p,e 

Re 

[15] 


Patel et al . (ref. 15) give a good description of these functions. 

The three functions given by Lam and Bremhorst (ref. 11) are as follows: 
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[18] 


This turbulence model has been chosen for use in this paper for a number 
of reasons. First, it has been used, by Rodi and Scheuerer (ref. 16), to 
predict transitional flow over a turbine blade with a boundary layer code. 
Second, when compared to other LRN k-e models it has been shown by Patel et 
al. (ref. 15) to be amongst the best in predicting the characteristics of 
fully turbulent flow. Third it has very simple boundary conditions and in 
general the model is clean and straightforward to code. 

Transformation of the Basic Equations 


The continuity equation [1] can be expanded to give, 


dijm) diftv) dip w) 
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3x 


dy 


d z 


[19] 
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The set of conservation equations describing the transport of momentum [2] 
can be written in Cartesian coordinates for a scalar variable as, 


d(p\xf>) dipv<j>) d(pv/(j)) _ d_ 


dx 


dy 


dz 


dx 


r <f>d <£' 

1 dx 


dy 


- dy. 


dz 




+ S 


<f> 


[ 20 ] 


where F ^ is an effective diffusion coefficient and denotes the source term. 
For the momentum equations the source term is a function of the pressure 
gradient. When new independent variables £,r) and f are introduced, Eq. [19] 
and [20] changes according to the following transformation £=£(x,y,z), 
r}=rj{x,y ,z) and £=f(x,y,z) as represented in Fig. 1. 


The transformation of partial derivatives is presented in the Appendix. 
Upon substitution and after some manipulation the basic equations [19] and 
[20] are transformed into the following form. 


d^ G ) dip G ) dipG ) 
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Vector Form of the Basic Equations 


Only the source term takes on a different form for the three different 
directions. These are, 

(1) <j> = n; S* = -£ 

(2) 4> = v; ^ - ^ 

(3) <f> = w; 

The source terms for the scalar equations considered such as the energy, 
turbulent kinetic energy and dissipation equations are mainly of scalar form. 

Discretization of the Transport Equations 
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The conservative form of the equation is derived for numerical 
computation, the integral form of the continuity equations [21] over the 
control volume, as shown in Fig. 2, is obtained as follows, 
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[24] 
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and the transport equation [22] is, 
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w 


where n, s, e, w, u and d are the locations of the intersection between the 
control volume faces and the grid lines. 


Equations [24] and [25] can be rearranged as follows: 

, d 


and. 


(pG^Af)® + (pG^Af)* + (pG^A^ = 0 
(J + (J/fif)" + (J { AfA,)^ = 

{^i»M ♦ ■$)**}! 

+ *MY^) 

+ H 3 ||) + S^A^A^AC 


[26] 


[ 27 ] 


where , 
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[ 28 ] 



which are the flux vectors through the faces of the control volume in Fig. 2. 

The source term is assumed to be acting at the centre, and constant within 
the control volume. 

The discretized form of the continuity equation [26] has been described in 
detail by Hobson and Lakshminarayana (ref. 10) for two-dimensional flows. 
Equation [27] is discretized by using central differencing and linear 
interpolation in the computational domain. The coefficients of the resulting 
difference equation are modified according to the hybrid scheme of Spalding 
(ref. 17). For a cell Reynolds number less than or equal to two the hybrid 
scheme uses central differencing including both the convection and diffusion 
terms. As the cell Reynolds number exceeds two the hybrid scheme switches to 
upwind differencing with only the convective contribution; the diffusion terms 
are neglected. 

Boundary Conditions 


No-slip boundary conditions are used for the velocities at a solid 
wall, where they are all set equal to zero. For cascade flows, 
periodic boundary conditions are implemented upstream and downstream 
of the blades on the periodically generated meshes. When solving 
along a periodic line in the blade-to-blade plane, two consequences 
result, which alter the tri-diagonal implicit matrix system of size 
N, which is usually used to solve between solid boundaries. First, an 
additional grid point is incorporated into the system, resulting in 
an N+l system of simultaneous equations. Second, a cyclic tri- 
diagonal matrix system results, which is solved with the algorithm 
presented by Napolitano (ref. 18). 

The inlet flow velocities are specified and fully developed flow is 
assumed to exist at the exit plane. The exit plane has to be situated 
far enough downstream of the blade row, where the assumption that the 
streamwise gradients are small is valid. 

For the pressure substitution method the boundary condition required to 
solve the pressure equation is that the normal derivative of the pressure 
vanishes at solid boundaries, which is a reasonable assumption for high 
Reynolds number. 

The boundary condition required to solve the pressure correction equation 
is that the normal derivative of the pressure correction, p’ , vanishes at the 
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boundary in the computational domain. The actual pressure value at the 
boundary is determined by extrapolation from the interior nodes to the wall 
with the condition of zero normal derivative of the pressure. 

Since Neumann boundary conditions are used for either the pressure 
correction method or the pressure substitution method, the pressure must be 
kept constant at one grid point. This control point is usually taken to be one 
of the inlet grid points, where the pressure gradients are small. 

For the solution of the energy equation to predict the variation of heat 
transfer coefficients over turbine blades, the blade wall temperature is 
specified from experimental results. 

The inlet turbulent kinetic energy and dissipation rate is determined by 
specifying an inlet turbulence intensity, Tu, and length scale, , associated 

with the inlet flow. Most sets of experimental data have the turbulence 
intensity quoted, however the length scale has to be determined from the 
physical dimensions of the turbulence generator. In most cascade tests, grids 
are used upstream of the cascade, to generate the required turbulence levels. 
The length scale is of the order of the characteristic dimension of the rods 
that make up the grid. The inlet turbulent kinetic energy is thus determined 
from the following, 


k i = f <*“ V* 


129] 


where tb is the inlet total velocity. The inlet dissipation rate is then 
determined from the following, 


£ . 


1 



[30] 


The boundary condition at the wall for the turbulent kinetic energy, k, 
takes on two definitions. For the "standard" form of the equations, with wall 
functions, the following condition is applied. 


5k 

^ wall 


[31] 


The value of k is set equal to zero at the wall in the LRN model. 


The idea behind the wall function method is that the turbulent length 
scale increases linearly, with distance from the wall beyond the viscous 
sublayer. The fluxes of momentum to the wall are then supposed, if the first 
grid point is in the fully turbulent region, to obey the relation, 


A-TC 1 / 4 k I/2 
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Here u , t , and y are, respectively, the time-averaged velocity of the 

fluid at the first grid point P away from the wall, the shear stress on the 

wall in the direction of the velocity u and the distance of the point P from 

P 


the wall. E is a constant and is a function of the wall roughness, 
approximately equal to 9.0 for a smooth wall, and /c (=0.42) is the von Karman 
constant . 


When calculating k , it is necessary to assign a value for the average 

energy-dissipation rate over the control volume, this is to be deduced from 
the assumption that, 




k 3/2 
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1 / 2 -, 


[33] 


The appearance of this logarithmic function results from the necessity to 
presume that, 


£ 

P 


(C l/2 k ) 3/2 
M P 


[34] 


which is the boundary condition used for the dissipation rate of turbulent 
kinetic energy. 


For the LRN model the boundary condition for the dissipation equation is, 



[35] 


It is clear that the boundary conditions for the low-Reynolds number model 
are simpler than wall functions, particularly in their implementation into a 
computer code. Because of the requirement that the viscous sublayer be 
adequately resolved, the LRN was used only for two-dimensional flows. Wall 
functions were used for the three-dimensional flow test cases as a result of 
limited computer memory. 

Solution Algorithm and Relaxation Technique 


The pressure substitution method (PSM) is used for the solution of two- 
dimensional flows and is fully described in reference 10. Upstream and 
downstream of the cascade the periodic solver according to Napolitano (ref. 

20) is used and between the blade surfaces a tri-diagonal matrix solver is 
used for the successive line under-relaxation (SLUR) of the computational 
domain. Implicit solution is in the cross-stream direction and the line sweeps 
are performed in the streamwise direction. 

The three-dimensional calculations were carried out with the pressure 
correction method. This algorithm is fully described by Rhie (ref. 12) and is 


292 



a logical extension of the two-dimensional form as described in reference 7. 
The alternating direction implicit (ADI) method with the implicit periodic 
solver is used to calculate the three-dimensional flow through an end wall 
cascade . 

The alternating direction implicit (ADI) method will be described by 
refering to Fig. 3. SLUR proceeds by taking all lines in the same direction in 
a repetitive manner. The convergence rate can be improved by following the 
sequence by rows, say, by a second sequence in the column direction. Such 
sequencing was used in the cross-stream plane for the cascade flows and this 
planar solution was advanced in the streamwise direction. 

If there are 1=1, NI grid points in the blade-to-blade direction and J=1,NJ 
grid points in the spanwise direction, as shown in Fig. 3, then on a cross- 
stream plane between the blades the ADI solution is performed from grid points 
2 to NI-1 and 2 to NJ-1 . However upstream (from grid points K=2 to KSTART) and 
downstream (from K=KFINIT to NK-1) of the cascade, as mentioned earlier, an 
additional periodic grid point is included in the blade-to-blade surface. Thus 
for the ADI line sweep in the blade-to-blade direction the periodic solver is 
used between grid points 1 and NI-1 (the values at NI then being set equal to 
those at 1). During the second ADI sweep in the spanwise direction, an 
additional line relaxation is performed on the periodic line, at 1=1, from J=2 
to NJ. Once again, to ensure periodicity the values at I=NI along the line 
J=2,NJ-1 are set equal to those at 1=1. This completes an ADI sweep in the 
cross plane of the computational domain representing a cascade geometry. The 
solution is now advanced in the streamwise, or K, direction. 

As the periodic solver is an extension of the Thomas algorithm (ref. 19), 
this solution procedure is believed to be consistent. The experience gained 
with the periodic solver in the solution of two-dimensional flows, showed that 
periodic solutions are efficiently achieved in an implicit manner, and thus 
was readily extended to three-dimensional flows. 

RESULTS AND DISCUSSION 


Laminar Flow in a Cascade 


The test case is that of flow through a cascade of symmetrical NACA 0012 
airfoils, which Rosenfeld and Wolfshtein (ref. 20) also computed. A comparison 
between the convergence histories of the pressure correction method (PCM) and 
the new pressure substitution method (PSM) was conducted for this 
conf iguration. 

The cascade had a solidity, or pitch-to-chord ratio, of 1.0, and the 
airfoils had a stagger angle of 30 degrees. The Reynolds number, based on 
chord length, of the inlet total velocity was 1000. The incidence angle was 
varied between -10 and 15 degrees. They used 81x37 grid points with a 
computational "0" type grid. Their formulation was based on the vorticity and 
stream function equations. 

The present computations were performed on a "H" type grid with 57x51 grid 
points. Almost the same number of grid points were used for the two studies, 
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however the "0" grid used by Rosenfeld and Wolfshtein (ref. 20) gives more 
grid points on the airfoil surface that the present "H" grid. The comparison 
between the PSM and PCM algorithms were performed at various relaxation 
parameters for the momentum equations. 

Fig. 4 shows the comparison of the convergence rate for the highest 
relaxation parameters (^ u = w y = 0.8, and u = 1.0). The PSM shows a significant 

improvement over the PCM which seems to indicate that at least the initial 
stages of convergence of the PSM is more favorably affected by the relaxation 
parameter than the PCM. The above test case was compared with zero degrees 
inlet incidence angle. A comparison of the predicted static pressure 
coefficient over the airfoil with that computed by Rosenfeld and Wolfshtein 
(ref. 20) is presented in Fig. 5. The absolute levels of pressure coefficient 
do not compare well. However, the predicted negative lift coefficient (C^= 

-0.0310) is within 4% of that predicted by Rosenfeld and Wolfshtein (ref. 20) 
who predicted C^= 0.0322. The negative lift force created at this incidence 

angle causes a negative turning angle of -0.7 degrees through the cascade. The 
predictions using the PSM and PCM resulted in plots on the same curve. 

An additional test case was computed with the PSM to predict the flow at 
the -10 degrees incidence case. The comparison with the predicted pressure 
coefficient as computed by Rosenfeld and Wolfshtein (ref. 20) is shown in Fig. 
6. This comparison is better than the zero degree incidence test case, and 
the predicted C^= -0.389 compares within 3% of their prediction of C^= 

-0.39711. The present method does not suffer from pressure field oscillations 
around the trailing edge. 


Turbulent Flow Through a Cascade of Double Circular Arc Profiles with Flow 
Detachment 


For the test case of turbulent flow through a cascade, the predictions 
made by the present method were compared with the experimental data obtained 
by Zierke and Deutsch (ref. 21). The blade section was a double circular arc, 
The computed test case was for the near design case of -1.5 degrees incidence. 

The modified version of the GRAPE code was used to generate the H-grid 
which extended half a chord upstream and one and a half chords downstream of 
the blade. The inlet angle of the grid was aligned with the incoming flow at 
51.5 degrees and the outlet grid angle was set equal to the measured outlet 
flow angle of 2.1 degrees. The Navier-Stokes calculations were on a 130 
streamwise by 100 tangential computational grid as shown in Fig. 7. The 
residuals were reduced by four orders of magnitude in roughly 2000 iterations. 
This corresponded to 20 minutes on the ETA-10 supercomputer, only using the 
optimization capabilities of the compiler which vectorizes inner DO-loops. The 
slope of the convergence plots were not monotonic and flattened out after 300 
iterations . 

Fig. 8 shows a comparison of the calculated and measured static pressure 
distribution. The experimental measurements indicated a possible separation 
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region near the trailing edge. Surface flow visualisation tests using the 
chemical sublimation method showed, with a 95% confidence level, a region of 
low shear stress at 45.1% ±2.3% chord. 

The comparison of the blade surface skin friction distribution is 
presented in Fig. 9. Transition of the pressure surface boundary layer from 
laminar to turbulent flow is predicted by the technique. The onset of 
transition is predicted to be at 10% chord, which is early compared to the 
experimentally determined transition process. Transition seems to be complete 
by 30% chord, thus the length of the transitional region is in agreement with 
that determined experimentally. The experimental points were determined by 
fitting a spline through the measured boundary layer profiles. The level of 
skin friction coefficient in the fully turbulent region is accurately 
predicted, as well as the final increase at the trailing edge. Good agreement 
is achieved between the predicted and experimental skin friction coefficient 
on the suction surface of the blade. Shown on this figure is the separation 
point as determined by the flow visualisation technique and the separation 
point as predicted by the code with good agreement. The levels of skin 
friction are very close to zero over most of the rear part of the blade. 

Turbulent Flow Through a Cascade of Turbine Nozzles Including Heat Transfer 


The computer program was used to compute turbine cascade flows for which 
Turner (ref. 22) conducted heat transfer measurments. A review of the 
literature showed that few incompressible flow experiments, including heat 
transfer measurements, have been documented. Turner (ref. 22) tested a nozzle 
blade section, as shown in Fig. 10, that was the same as the one used by 
Bayley and Turner (ref. 23) and Bayley and Wood (ref. 24). The blade chord was 
70 mm and the aspect ratio and pi tch-to-chord ratio were 1.34 and 0.65, 
respectively. The stagger angle was 40.6 degrees, and the inlet and outlet 
flow angles were 0.0 and 61.0 degrees, respectively. Three exit Mach numbers 
(0.75, 0.65, 0.55) were considered. At each Mach number tests were done at 
three turbulence intensities, namely 5.9%, 2.2% and 0.45%. 

The cascade tunnel delivered 1.5 atmospheres of air, at 90 degrees C, to 
the test section. The blades were internally cooled with air at 20 degrees C, 
which produced a sufficient temperature gradient for heat transfer 
experiments. At the highest turbulence level the average gas-to-wall 
temperature difference was 28 degrees C and 40 degrees C at the lowest 
turbulence level. The exit Reynolds number based on the nozzle chord was 

1 .05x10? 

Since the present predictions are valid for incompressible flow, only the 
low Mach number test case was considered. For this test case, the aerodynamic 
experimental test data were obtained from Bayley and Wood (ref. 24). Although 
porous blade sections were tested, they did present blade surface velocity 
distributions for the solid blades. 

The turbulence of the flow at entry to the cascade was varied by placing 
various grids and drilled plates one chord distance upstream of the blade. 
Unfortunately, neither the diameter of the grids nor their character 
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dimensions were quoted by Turner (ref. 22). The turbulence length scale, , 

was thus set equal to the leading edge radius (0.26cm) of the nozzle blades 
for all the test cases considered. This value was used in equation [30] to 
determine the inlet levels of the dissipation rate e. 

As shown in Fig. 10, a 100x132 computational H-grid was used with the 
normal grid spacing at the blade surface such that the first grid point was at 

a y + of 1.0. The inlet angle of the grid to the cascade was set equal to zero 
degrees, and the exit grid angle was set equal to 60 degrees. 

The predicted velocity distribution over the blade surface is compared to 
the distribution calculated from the static pressure measurements in Fig. 11. 
The local velocity was derived by assuming isentropic expansion from the inlet 
total pressure to the static pressure, either measured or predicted, recorded 
on the blade surface. The sudden acceleration of the flow around the leading 
edge on the suction surface is captured, and the subsequent acceleration on 
the front section of the blade. The maximum velocity level is accurately 
predicted, and the diffusion process that results on the suction surface after 
the turbine throat is adequately resolved. The level of local velocity is 
slightly high on the rear section of the suction surface. The rapid diffusion 
process on the rear part of the suction surface results in flow separation 
which is predicted at the trailing edge. The continual flow acceleration 
process on the pressure surface is well predicted over its entire length. 

Turner (ref. 22) showed that the mean heat transfer results for a turbine 
blade operating at one Reynolds number, but at different turbulence levels, 

vary substantially. At a Reynolds number of l.OxlO 6 the mean heat transfer 

coefficient varied from 8.5x10 to 1.05x10 for varying inlet turbulence 
intensity levels from 0.45 to 5.9% respectively. For different turbulence 
levels the slope of the Nusselt number versus Reynolds number also varied, the 
actual slopes are presented by Turner. A distinct increase in the blade 
surface temperature was measured. No data were presented on the cooling system 
heat transfer, thus a complete energy balance could not be determined for the 
cascade. The measured distribution of surface temperature at 0.45 and 5.9%, as 
shown in Fig. 12, was used as temperature boundary conditions. Linear 
interpolation was used between these two levels to determine the temperature 
boundary conditions at 2.2% turbulence intensity. 

An example of the prediction of the variation of local heat transfer 
coefficient over the surface of the blade is presented in Fig. 13 for varying 
inlet turbulence intensities. As mentioned earlier the stagnation point is 
situated on the pressure surface, this is illustrated by the maximum heat 
transfer rates in this region which is offset toward the concave surface. The 
levels of heat transfer are in good agreement with the experimental data in 
the stagnation point region. No sudden transition from a purely laminar to a 
purely turbulent flow is observed, and this trend is well predicted by the 
program. At an inlet turbulence intensity level of 5.9% the predictions show 
some form of transition of the boundary. Although this is not in complete 
agreement with the experimental data, this could be due to the mixing 
augmentation on the concave surface dominating the favorable pressure 
gradient. It is well established that in a two-dimensional boundary layer, 
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flow over a convex surface damps mixing, while a concave surface leads to 
augmentation as described by Bradshaw (ref. 25). On the pressure surface the 
prediction of the increase in the average heat transfer coefficients for 
increasing turbulence levels is very good. This is surprising, considering 
that a constant turbulent Prandtl number (=0.9) was used for these 

calcuations. The average gas-to-wall temperature difference was 30°C; however, 

( T / T rvy 

the gas-to-wall temperature ratio gas wall ~ 1.09) was representative of 
those found in the uncooled turbine stages of a gas turbine engine. 

For the suction surface, a sharp increase in heat transfer between 80% and 
90% chord is apparent for the low turbulence levels, and the increase is 
between 50% and 90% for the high turbulence level case. Fig. 11 indicates 
that on this surface a favorable pressure gradient, decreasing in magnitude, 
extends 60% of the chord downstream of the leading edge. The boundary layer 
then experiences an adverse pressure gradient and the high heat transfer rates 
in this region are indicative of transition from laminar to turbulent flow. 

For the 0.45% and 2.2% test cases the transition process was experimentally 
determined to be abrupt, which the LRN turbulence model is not able to 
accurately simulate. Once again the model, although only for low turbulence 
levels, has tended to smear the transition process over a longer distance and 
the onset of transition is predicted to occur too close to the leading edge. 
The experimental data showed little or no effect for increasing turbulence 
levels from 0.45% to 2.2% over the front 70% chord of the suction surface. A 
further increase in turbulence to 5.9% produced a significant increase in heat 
transfer. The length of transition of the boundary layer for this case is 
longer than the low freestream turbulence intensity test cases. The transition 
process is adequately predicted for this high turbulence level case, which is 
the one of the strong points of the LRN k-e turbulence model. The prediction 
showed an increase in heat transfer from 0.45% to 2.2% which is not in 
agreement with the experimental data. This could be due to the interpolated 
blade surface temperature that was used as the boundary conditions for the 
2.2% turbulence case. The actual blade surface temperature, which was not 
documented, could be different from that used. This could explain the 
discrepancy between these results. 

Flow in a Turbine End Wall Cascade 


The final three-dimensional flow test case considered is secondary flow 
and losses in a turbine cascade. Gregory-Smi th and Graves (ref. 26) 
investigated the flow in a cascade of large scale rotor blades of some 110 
degrees of turning. The flow was traversed with 3-hole and 5-hole cobra type 
pressure probes. Hot-wire anemometry measurements were also carried out in the 
linear cascade consisting of blades scaled from the midspan section of a 
high-pressure turbine rotor design. The main geometric parameters for the 
cascade are given in Table 1. The blade profile coordinates were obtained from 
Gregory-Smi th (ref. 27). 
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Table 1 Cascade Geometric Parameters 


Inlet flow angle 

42.75 

deg. 

Inlet blade angle 

52.25 

deg. 

Exit flow angle 

-66.60 

deg. 

Exit blade angle 

-67.50 

deg. 

Blade chord 

216 

mm 

Blade axial chord 

175 

mm 

Span 

457 

mm 

Pitch 

191 

mm 

Reynolds number 

(based on chord and exit velocity) 

5 x 10 5 



The estimated inlet (99%) boundary layer thickness was 102.0 mm, the 
displacement thickness was 12.1 mm and the momentum thickness was 9.9 mm. 

These data were obtained at the first slot which was 14% of an axial chord 
upstream of the cascade. The measured inlet free-stream turbulence intensity 
was 1.4% which was due to the honeycomb flow straightener placed upstream of 
the cascade. The inlet boundary layer profile agrees well with a 0.12 power 
law profile. This profile was used for the input velocity distribution during 
the computational analysis. 

The three-dimensional "H" grid was generated from a two-dimensional 
blade-to-blade grid which was projected in the spanwise direction. Grid 
clustering was used at the endwalls, at the blade sufaces and around the 
leading edge and trailing edge of profile. In the cross-plane, 31 x 31 grid 
points were used with 81 in the streamwise direction. A three-dimensional view 
of the grid is presented in Figure 14. Although the problem is symmetrical 
about the mid-span plane, the computation was performed for the complete 
cascade to check the code’s ability to compute a symmetric solution. 

The three-dimensional flowfield was initialised with the solution of the 
two-dimensional flow through the cascade. The two-dimensional version of the 
code was run to 500 iterations, giving two orders of magnitude convergence. 

The velocity field showed most of the correct features of two-dimensional flow 
through a cascade. As the LRN turbulence model was used (on this coarse grid) 
the boundary layer growth is most probably too rapid; however, this was 
considered to be a good initial guess for the three-dimensional flow. 

The velocity field was scaled in the spanwise direction according to the 
inlet boundary layer profile at the endwalls. The pressure distribution was 
assumed constant in the spanwise direction. Less than one order of magnitude 
convergence was obtained for this case, which took 4000 CPU seconds on the 
ETA-10 at the John von Neumann Computer Center. Because of limited computation 
time the solution was not continued further. 

The main objective of this test case is to attempt to qualitatively 


298 




predict the three-dimensional viscous phenomena inside a turbine blade 
passage . 

Fig. 15 shows the velocity vector plot for the end wall region. Flow 
reversal occurs at the leading edge which is associated with the roll-up of 
the inlet boundary layer at the leading into the so-called horse-shoe vortex. 
The subsequent movement of this vortex is not captured as the secondary flow 
convects the pressure side vortex over to the suction surface of the adjacent 
blade too rapidly. 

The horse-shoe vortex and passage vortex begin in the leading region close 
to the endwall. An "H" grid does not allow the adequate resolution of the 
severe flow gradients in this region. This region should be resolved with a 
"C" grid which is continuous around the leading edge. A higher order 
turbulence model needs to be considered if the finer details of this complex 
flow is to be adequately resolved. 

CONCLUSIONS 


Various test cases have been presented. These show that in two dimensions 
the pressure substitution method is superior to the pressure correction 
method. This method is ideally suited for the computation of incompressible 
flows . 

For complex flow situations such as turbulent flow through a compressor or 
turbine cascade, the code exhibits satisfactory convergence behavior. The 
global parameters such as pressure distribution, lift coefficient and heat 
transfer are well predicted within engineering accuracy. The agreement between 
the computed and measured blade surface skin friction coefficients and heat 
transfer coefficients is very good. Details of the leading edge flow and 
particularly the trailing edge separation have been resolved, which shed more 
light on these complex flow regions. 

For flow over curved surfaces with large pressure gradients the minimum 
turbulence model that must be used is the low-Reynolds number k-£ model. Its 
ability to predict transition dependent on freestream turbulence intensity is 
well known, and it was used in the present study to predict separated flow and 
flow through a turbine nozzle including heat transfer. 

Experimental investigations of turbulent flow through cascades should 
include detailed measurements of the inlet freestream turbulence intensity and 
turbulence length scales. These two parameters are of vital importance for the 
correct specification of the inlet boundary conditions. 

The three-dimensional turbulent flow through a turbine end wall cascade 
has been qualitatively resolved. To increase the accuracy of the predictions 
it is recommended that a higher order turbulence model be used which can 
simulate the laminar flow leading edge region. 
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APPENDIX A 


Partial derivatives of any function are transformed according to, 

§_ _ * d_ d_ ^ 

dx ^x d£ + ^x dr] + ^x dg 

d_ _ t d_ <L , JL 

dy S + ^y dr) + ^y dg 

d t d_ d_ d 

dz ^z d£ + \ dr) + 

where , 
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and , 


£ = (y z. - y .z )/I 

x VJ r? f r? 

{ y ^‘Vf ~ hV' 1 
fz = <Vf “ W" 


"x ' V{ )/Ir 

*7 y = (x £ z f “ Vf )/J 

\ =- (x e y ? - V( )n 

7 = (y f 2 >, - v« )/J 

C =-(x>.z - x z„)/J 

b y £7 v £ 

C = (x*y - x y*)/J 
5 z ^77 77 £ 


J = 


x £ x >? x <r 
y e y 7 y c 
z £ z v z f 


which is the Jacobian of the transformation. 

The contravariant velocities are defined as, 

G = £u + £v + £w 

1 s x s y s z 

G = 7?U + 77V + 7?W 

2 'x 'y 'z 

G„=Cu+Cv+fw 

3 3 x 3 y 3 z 

and the additional transformation metrics as, 

>2 >2 >2 

a = ( + f + £ 

^x ^y s z 

a 2 2 2 

^ = \ * ”y + "z 

7 = f 2 + f 2 + f 2 

' 3 x 3 y 3 z 


and, 


H, = ^ 77 + £ 77 + £ ?7 

1 s x'x 77 s z'z 

H = £ s + £ C + £ f 

2 s x 3 x 77 s z 3 z 

H = 77 C + )) f + W C 

3 'x 3 x 'y 3 y 'z 3 z 
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Fig. 1 Transformation Representation 


Computational Domain 

Fig. 2 Control Volume 



Fig. 3 Alternating Direction Implicit Relaxation 


of a Three-Dimensional Cascade Geometry 


303 




ITERATION NUMBER 

Fig. 4 Comparison of Convergence Rates for 

PSM and PCM; Laminar Flow in a Cascade, 

w =w = 0,8, w =1.0 
u v ’ p 



I'ig. 5 Comparison of Predicted Pressure Coefficient 

for a Straight Cascade of NACA 0012 Blade Profiles 
at 30 Deg. Stagger and 0.0 Deg. Incidence 



Fig. 6 Comparison of Predicted Pressure Coefficient 
for a Straight Cascade of NACA 0012 Blade Pr 
at 30 Deg. Stagger and -10.0 Deg. Incidence 
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Fig. 10 Computational Grid for the Turbine Nozzle Cascade 
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Fig. 11 Comparison of the Predicted and Measured Local Velocity 
Distribution for the Turbine Nozzle 
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Fig. 12 Measured Blade Surface Temperature Distribution: Turner (ref. 22) 
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Fig. 14 Computational Grid for the Turbine End Wall Cascade 
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Fig. 15 End-Wall Velocity Vector Plot Showing the Development of the 
Horse Vortex System 




